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ABSTRACT 

Although the long-term numerical integrations of planetary orbits indicate that our 
planetary system is dynamically stable at least ±4 Gyr, the dynamics of our Solar 
System includes both chaotic and stable motions: the large planets exhibit remarkable 
stability on gigayear timescales, while the subsystem of the terrestrial planets is weekly 
chaotic with a maximum Lyapunov exponent reaching the value of 1/5 Myr -1 . In this 
paper the dynamics of the Sun- Venus-Earth-Mars- Jupiter-Saturn model is studied, 
where the mass of Earth was magnified via a mass factor ke- The resulting systems 
dominated by a massive Earth may serve also as models for exoplanetary systems that 
are similar to our one. This work is a continuation of our previous study, where the 
same model was used and the masses of the inner planets were uniformly magnified. 
That model was found to be substantially stable against the mass growth. Our simu- 
lations were undertaken for more then 100 different values of ke for a time of 20, in 
some cases for 100 Myrs. A major result was the appearance of an instability window 
at ke ~ 5, where Mars escaped. This new result has important implications for the 
theories of the planetary system formation process and mechanism. It is shown that 
with increasing ke the system splits into two, well separated subsystems: one consists 
of the inner, the other one consists of the outer planets. According to the results the 
model became more stable as ke increases and only when ke ^ 540 Mars escaped, 
on a Myr timescale. We found an interesting protection mechanism for Venus. These 
results give insights also to the stability of the habitable zone of exoplanetary systems, 
which harbour planets with relatively small eccentricities and inclinations. 
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1 INTRODUCTION 

The determination of the stability of our Solar System is one 
of the oldest problems in astronomy. The question has been 
debated over more than 300 years, and has attracted the 
attention of many famous mathematicians over the course 
of history. The problem played a central role in the devel- 
opment of non-linear dynamics and chaos theory. Despite 
the considerable efforts, we do not possess a definite an- 
swer to the question of whether our Solar System is sta- 
ble or not. This is partly a result of the fact that the def- 
inition of the term stability is not unambiguous when it 
is used in relation to the problem of planetary motion. In 
addition to the vagueness of the concept of stability, the 
planets in our planetary system show a character typical 
of dynamical chaos. The physical basis of this chaotic be- 
haviour is now partly understood as a consequence of reso- 
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nance overlapping and three body resonances fl Lecar et alJ 
l200ltlMurrav fc Holmanll999l:lNe"svornv fc Morbidellll999¥ 

which can manifest themselves in dramatic and relatively 
sudden changes in an orbit. In the last two decades sev- 
eral numerical stability studies were performed in order to 
throw light on the question. At pres ent the longest numerica l 
integrations published are those of llto fc Tanikawal i2002f) . 
where six long-term numerical integrations of all nine plan- 
ets, covering a time-span of several 10 9 and 10 11 years are 
discussed. Their fundamental conclusion is that the Solar 
System seems to be stable in terms of the Hill-criteria at 
least over a time-span of ±4 Gyr. Moreover it turned out 
that during the integration period all the planetary orbital 
elements have been confined in a narrow region. 

On the other hand accord ing to Laskar' s semi-analytical 
secular perturbation theory jLaskarlll98ll . the terrestrial 
planets', especially Mercury's and Mars' eccentricities and 
inclinations show large a nd irregular variations on a time- 
scale of several 10 year (lLaskarlll996l) . 
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Figure 1. Eccentricities vs. the semimajor axes of the observed 
cxtrasolar planets. The x-axis is logarithmic. The positions of the 
Earth, Jupiter and Saturn are also indicated as diamonds with a 
plus sign in the middle. 



Nowadays to study the stability of the Solar System, or 
its variants, as a representative of the different planetary sys- 
tems has become part of the frontline research. Over the past 
few years the detection of planets outside the Solar System, 
the so called exoplanets, has greatly stimulated the stabil- 
ity studies of planetary systems. New exoplanets are being 
discovered on a regular basis; more than 150 (April, 2005) 
exoplanets are now known. There are 136 systems consisting 
of a central star and a gaseous giant planet, and 14 multiple 
systems with two, three and four planets. The so far discov- 
ered exoplanets have a minimum mass range (m • sin(i p )) 
from 0.042 m.j to 17.5 mj where mj is Jupiter's mass and 
i p is the inclination of the orbital plane with respect to the 
plane of the sky. Since i p is unknown a precise mass cannot 
be determined, only a lower mass limit. Because of the tech- 
nical limitations only planets of Neptune mass or above can 
be detected and then only if they are less than 5 AU or so 
from the star . Therefore more than 90 % of these planets 
are orbiting their host star well inside Jupiter's orbit. There 
are major differences between the characteristics of the so 
far observed systems and those of the Solar System. Most 
of the planets have minimum masses substantially greater 
than that of Jupiter - up to six or even more times the 
mass of Jupiter. Dozens of planets are orbiting very close to 
their hosting star, with semimajor axes down to 0.04 AU. 
Finally, planetary orbits are found with large eccentricities, 
up to approximately 0.7, plus a few greater, significantly 
greater than the highest eccentricities observed for planets 
in our Solar System. These characteristics are depicted in 
Fig. where the planets' eccentricities are plotted against 
their semimajor axes, the locations of the Earth, Jupiter and 
Saturn are marked as diamonds. From Fig. Q it is apparent 
that our planetary system may serve as model case for those 
planetary systems which have small eccentricities. Presum- 
ably these are also the ones where we may expect stable 
terre strial planets moving in habitable zones iAsehari et alJ 
l2004h . 

In previous papers JPvorak fc Sulill2002l; iDvorak et ail 

120051) the dynamical evolution of a simplified Solar system 
was studied. The model consisted of the Sun, the three most 



1 With OGLE, which is based on optical gravitational lensing it 
is possible to detect exoplanets with only Earth-mass. 



massive terrestrial planets (Venus, Earth, Mars), Jupiter 
and Saturn. The masses of the inner planets were uniformly 
magnified by a mass factor n. It turned out that the dif- 
ferent systems remained stable up to 10 Myr for k ^ 220. 
Stable islands were found for k = 245 and 250, which is 
a well-known property in such regions which are close to 
the last stable orbit in the chaotic domain. We have shown 
that the dynamical coupling of Venus and Earth and that 
of Jupiter and Saturn remained unbroken for all studied k. 
On the other hand the motion of Mars was not coupled to 
any other planet, what may be a reason for the fact that 
it was always Mars which caused the decay of the system 
after close approaches with Earth. However the remarkable 
stability of these model planetary systems suggests that ex- 
oplanetary systems with configuration like our Solar System 
may harbour moderately or even very massive terrestrial-like 
planets. 

In the present work our aim was to study and ana- 
lyze the dynamical evolution and the stability of the system 
with respect to the masses involved. Contrary to our previ- 
ous study, in the present work only the mass of the Earth 
was increased via a mass factor Kb- Furthermore, the ex- 
amined systems may be considered as models for individual 
exoplanetary system, and the results can be applied to them. 
Section 2 explains our dynamical model, the applied meth- 
ods and section 3 is devoted to a detailed description of the 
results in the ke G [1,600] region. Finally, we discuss the 
results and the implication for exoplanetary systems. 



2 DESCRIPTION OF THE MODEL AND 
METHODS 

The applied dynamical model consisted of the Sun, Venus, 
Earth, Mars, Jupiter and Saturn. We have chosen this model 
for two reasons. To speed up the numerical integrations, we 
have omitted Mercury, Uranus and Neptune. Because of the 
small mass of Mercury, it only slightly perturbs the motion 
of the terrestrial planets and it has a short orbital period 
which would require a reduction in the integration step-size 
resulting in increased CPU time. Even though Uranus and 
Neptune are massive planets, they are evolving around the 
Sun more than twice as far as Saturn, and so they do not in- 
fluence the motion of the inner planets significantly. These 
simplifications are also justified by the fact, that the dy- 
namics of the inner planets are dominated by Jupiter and 
Saturn. Furthermore the so far observed exosystems har- 
bour at most three of four planets. The above modification 
of the Solar System gives grounds for the adjective "simpli- 
fied" . Hereafter we will refer to our simplified Solar System 
as S 3 . The other reason was to make the model more like 
an exoplanetary system. The mass of the Earth was magni- 
fied via a mass factor ke (the masses of the other planets 
were unaltered), which resulted in such hypothetical plan- 
etary systems, whose characteristics agree very well with 
several exosystems (for example if we multiply the semi- 
major axes of the planets of 47 Ursae Majoris by 2.5, the 
resulting configuration is something like the Sun, Jupiter 
and Saturn system, if k_b = 1). According to these modifica- 
tions, our models are parallel with those exosystems, which 
harbour planets with small eccentricities and inclinations. 
In Fig. the masses of the planets are plotted as a func- 
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tion of ke- At he ~ 90 Earth is as massive as Saturn and 
when he ~ 300 Earth is as massive as Jupiter. The mass 
of the other planets are plotted by horizontal dashed lines. 
The mass distribution of the so far discovered exoplanetary 
systems is also depicted in Fig. [5] up to mass parameter = 
2. The integrations were done in the ke £ [1,600] region, 
which approximately corresponds to < mass parameter ~ 
2 interval, and therefore the mass distribution in Fig. [5] is 
displayed only in this interval. 

We have considered the S 3 as a non-linear Hamiltonian 
system, governed only by classical Newtonian gravitational 
forces between the objects of the model. The planets have 
been taken as point masses, and the Earth's Moon was not 
included in the models. The initial planetary orbital ele- 
ments and the actual masses are listed in Table Q 

To check whether the different S 3 setups belonging to 
different kes could be stable over a long time interval we 
utilized the very precise numerical integration scheme, the 
Lie-integrator method. This method is based on the inte- 
gration of differential equations with Lie-series and uses the 
property of recurrence formulae for the Lie-terms. The de- 
tails o f the method are describ ed in iHanslmeier fc Dvorak! 
il984l) and iLichteneggeil Jl984T> . The scheme is particularly 
effective in the case of highly eccentric orbits. The accu- 
racy of this integration technique is based on an automatic 
stepsize control, and it has been checked in several com- 
parative test computations with other integrators. Although 
symplectic integrators are very effective when eccentricities 
remain small, but the Lie-integrator is a better choice in 
studies, where very large eccentricity orbits may occur. 

The length of the integrations was fixed at 20 Myrs, in 
some cases at 100 Myrs, which was a trade-off between too 
long CPU time and the quality of the results. Since our inter- 
est focuses primarily on the inner three planets' motions, for 
which the orbital time-scales are much shorter than those of 
the outer two giant planets, the 20 Myrs seems a justifiable 
choice, alth o ugh it is known from earlier work (see for e.g. 
IJones et all (1200 J) ; Hones fc Sleet] te002l) 1 that exosystems 
can stay stable for hundreds of millions of years and then 
fall apart. In the case of weak chaos therefore the 100 and 
especially the 20 Myr timespan might be short and in some 
of the integrations where the terrestrial planets survived for 
20 or 100 Myrs, they might not survive significantly longer. 
For the sake of a comparison all results were derived on 
the same computer. In some cases we have performed com- 
parative integrations on different platforms (to obtain more 
information of the particular system). 

The conservative definition of the point at which sys- 
tems become unstable is when close encounter between two 
planets happens: two bodies approach one another within 
an area of the larger Hill radius. The consequence of such 
an event is the dramatic changes in the orbital elements of 
the two planets, and usually the escape of the planet with 
smaller mass. In this paper we define a system unstable, 
when an orbit crossing or a close encounter happens. This 
definition is somewhat more general and the instability can 
be directly connected to the eccentricity via the perihelia 
and aphelia distances. Henceforward we state that our model 
is dynamically unstable if orbit crossing or close encounters 
happen in the course of the integration. Using both criteria 
is clearly an extension of the conservative definition of in- 
stability. It is justifiable to incorporate the orbital crossing 



criteria in the definition since we know from experience that 
an orbital crossing in general leads to a close encounter in 
very short time. The main difference between the two defi- 
nitions is the time-scale of instability. We note that orbital 
crossing does not lead to close encounter in all cases when 
certain resonances have adjusted the planetary motions in 
such a way that the planets avoid each other. For example 
this is the case in the Neptune-Pluto pair. In our models 
no such protection resonances are present, henceforward we 
use the above definition to distinguish between stable and 
unstable S 3 setups. The integrations were not stopped after 
one of the above criterion had been met, but were continued 
until the integration timespan was reached. 

For an indication of stability we used a straightforward 
check based on the eccentricity. This osculating orbital el- 
ement shows the probability of orbital crossing and close 
encounter of two planets, and therefore its value provides 
information on the stability of the orbit. We examined the 
behavior of the eccentricities of the planets along the inte- 
gration, and used the largest value as a stability indicator; 
in the following we call it the maximum eccentricity method 
(hereafter MEM). This is a reliable indicator of chaos, be- 
cause the overlap of two or more resonances induce chaos 
and large excursions in the eccentricity. We know from ex- 
perience, that instability comes from a chaotic growth of 
the eccentricity. This simple check has already been used in 
other stability studies, and was found to be qu ite a powerful 
indicator of the stability c haracter of an orbit iDvorak et alJ 
l2003tlAsghari et allliooi) . 

2.1 The Laplace-Lagrange secular theory 

In order to find a theoretical explanation for the decay of 
the system we have applied the Laplace-Lagrange first or- 
der secular theory. This linear theory yields accurate results 
under the following assumptions: 

(i) no mean motion commensurabilities, 

(ii) no orbit-crossing, and 

(iii) the eccentricities and inclinations are small enough. 

Our models meet these criteria. However a complication 
is that the orbits of Jupiter and Saturn are close to a 
5:2 commensurability. Since the appearance of the Laplace- 
Lagrange theory, as a first approximation it has been exten- 
sively used in the studies of motion of the plane ts and of 
other Solar System bodies. In several researches llKnezevid 
Il986t lLaskarlll988l) the results of the first and higher or- 
der secular theories were compared. According to these 
studies the secular frequencies calculated from the Laplace- 
Lagrange theory are sufficiently accurate for our present re- 
search goal. The main discrepancies are in the case of Jupiter 
and Saturn. As we are interested only in the motions of the 
inner planets, therefore we will apply the Laplace-Lagrange 
theory to our models. From the above assumptions it is ap- 
parent that the precision of the theory does not depend on 
the planetary masses, accordingly there is no theoretical lim- 
itation on the mass factor's magnitude. 

Since the eccentricities and inclinations may vanish at 
remote epochs, it is better to use the Lagrange orbital ele- 
ments: 

/ h \ _ siniz? / p \ _ . sinfi . . 

\ k J cos w ' I q J cos 
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Figure 2. The mass distribution of the known exoplanetary systems are plotted together 
with the mass of the Earth as a function of the mass factor (kb) (solid line). The masses of 
the other planets were left unchanged, and they are represented by straight, dashed lines. The 
lower x-axis is the mass factor, the upper x-axis is the mass parameter of the exoplanetary 
system. The left y-axis is the mass of the planets in a logarithmic scale, the right y-axis is 
the number of planets. The width of the bin is 0.2 (this refers to the top axis) 



Table 1. Planetary orbital elements (JD 2449200.5) with respect to the mean ecliptic and equinox J2000. 
The quantities a, e. i, u>, Q and M denote the semimajor axis, eccentricity, inclination, argument of perihelion, 
longitude of ascending node and mean anomaly. In the last row the masses of the planets in Solar mass units 
are listed (Msm = !)• 





Venus 


Earth 


Mars 


Jupiter 


Saturn 


a 


0.723328 


0.999999 


1.523614 


5.202627 


9.545509 


e 


0.006747 


0.016716 


0.093443 


0.048370 


0.052420 


i 


3.394820 


0.000545 


1.850191 


1.304638 


2.485620 


OJ 


54.847892 


113.611521 


286.492727 


275.222227 


338.025839 


n 


76.691772 


349.288391 


49.573832 


100.470086 


113.651098 


M 


135.521541 


78.172620 


185.208769 


233.733076 


259.852365 


1/mi 


408 523.71 


332 946.047 


3 098 708.0 


1047.348 


3497.898 



Eq. (y associates the h, k, p and q Lagrangian-elements to 
the e, i, m, and Keplerian-elements, where e denote the 
eccentricity, i the inclination, w the longitude of perihelion 
and fi the longitude of ascending node. Using these vari- 
ables the general solution of the differential equations for 
the planets takes the following form: 

(£) - t^Z <"' + A>. ( 2 ) 
j=i 

3=1 

where the s index denotes the planet, the j index denotes 
the mode, N is the number of planets, and Lg are the 
amplitudes, gj and fj denote the secular frequencies and f3j 
and jj are the angular phases. 

The planet's orbital elements are described by Eq. 
and Eq. 10 , which are the sum of harmonic oscillations. Us- 
ing these formulae it can be calculated that the planet's ec- 
centricities and inclinations are varying between given limits 
with quasiperiodic oscillations. Due to the positive gj secular 
angular velocities the apsidal lines of the planets are rotat- 
ing directly, whereas the nodes accordingly to the negative fj 



secular angular velocities are rotating indirectly. Upon these 
mean rotations quasiperiodic variations are superimposed. 
Both of the apsidal and nodal motions can be approximated 
by average angular velocities, which are to a first approxima- 
tion equal with the frequencies of those harmonious terms 
which are multiplied by the largest amplitudes: 

sincc7 s . ,c ;i sin , _ , ... 

e s ■ « Afj J) (gjt + fa) (4) 

COSG7 s COS 

where M a (J) = max, |M S 0) | and the average angular velocity 
of the s'th planet is given by gj. 

In this manner on the basis of the amplitudes each sec- 
ular frequency can be associate with each planet. This asso- 
ciation is not mutually unambiguous. It may happen to as- 
sociate a certain secular frequency to more then one planet. 

The use of the linear theory seems to be in contradic- 
tion to the large eccentricities and inclinations that may be 
reached by a planet during the simulation. We emphasize 
that the forementioned theory was used at the beginning of 
the integration, when the inclinations and the eccentricities 
are small, and the above assumptions are therefore fullfilled. 

The above described linear secular theory was imple- 
mented in the MAPLE computer algebra program. With 
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Table 2. Summary of the integrations. 



ke interval 


Ak e 


Earth mass in 


Total time 






Jupiter mass 


[Myr] 


1 - 25 


1.0 


[1/300 - 1/12] 


20 


4-6 


0.1 


[1/75 - 1/50] 


100 


30 - 200 


5.0 


[1/10 - 2/3] 


20 


210 - 300 


10.0 


[2/3 - 1] 


20 


330 - 600 


30.0 


[1-2] 


20 



the aid of this application the formulae of the theory can 
be evaluated in a few seconds and it gives the complete first 
order solution of the problem. 



3 RESULTS 

In this section we give a detailed description and overview 
of the results of our simulations. The model was integrated 
for more then 100 different values of ke- The interval of 
data output was 100 years, the total amount of data is ap- 
proximately 6 GBs, and the total used CPU time is more 
than several thousands of days. All of the integrations were 
done on two Sun Fire 15000 supercomputers with 72 US- 
111+ 1200 MHz processor in each computer. In Table |2] the 
first column lists the different ke intervals, the second shows 
the stepsize (Akb) in the mass factor. A small stepsize is 
taken for 4 ^ he ^ 6 in order to explore the interesting 
behaviour of the system leading - very surprisingly - to the 
escape of Mars. In the third column we list approximately 
the mass of Earth in Jupiter's mass unit. The last column 
gives the time-span of numerical integrations. 

3.1 The 1 < n E «S 25 region 

As we expected, the orbital motions of the planets indicate 
long-term stability in most of our numerical experiments: 
no orbital crossings nor close encounters between any pair 
of planets took place in the course of the integrations. How- 
ever, a suprising result was the discovery of the instability 
window for the ke £ [4, 6] interval: at several ke values Mars 
escaped; the details can be found in the next subsection. 

Fig. depicts the results of the MEM. The maximum 
eccentricities (hereafter ME) of Jupiter and Saturn are ac- 
tually constant, with a value of 0.06, and 0.088, respec- 
tively (they are not shown in Fig. |3J. As one can clearly 
see from Fig. |3] the MEs of the Earth and Venus are rela- 
tively small, and both curves show a similar behaviour as a 
consequence of the well-known coupling between them. We 
note, that the ME of Earth for ke = 1 is 60% greater than 
for ke ~ 2. After some oscillations of Venus' and Earth's 
MEs, they stay almost constant with a value of 0.041 and 
0.035, respectively. In turn, the ME of Mars steadily in- 
creases with ke, and at ke = 5 it suddenly reaches a very 
high value, eM ars = 0.26091 (the perihelion distance of Mars 
is qMars = 1.126)! After this peak Mars' ME drops down to 
its starting value 0.1231, and begins to increase slowly and 
gradually with ke- 

The character of the variation of the planetary orbital 
elements does not change substantially over the course of 
the simulation, except for ke around 5. The variations of 
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Figure 3. Results of the MEM in the 1 < k e < 25 region 
(Ake = !)• The x-axis is the ke, the y-axis is the maximum 
eccentricity (ME). 

the semimajor axes are in the order of 10 -4 -10 -3 AU for 
the inner planets, and for Jupiter and Saturn they are ap- 
proximately 10~ 3 , 10~ 2 AU, respectively. The variations in 
Jupiter's and Saturn's semimajor axes are one order of mag- 
nitude larger than those of the inner planets, which is a 
consequence of the 5:2 near commensurability between their 
motion. For all simulations in this mass factor region the 
behaviour of the semimajor axes did not change, and their 
excursions remained in the above range. 

The eccentricities as a function of time of the inner plan- 
ets for ke = 1, 2, and 15 are shown in Fig. [I] (for the sake 
of a comparison the eccentricities for the actual masses are 
also depicted). As ke goes from 1 to 2, there is a substantial 
change in e for all three terrestrial planets. According to the 
MEM analysis (see Fig. EH , the MEs of Earth and Venus de- 
crease, which manifests themselves in greater perihelia and 
smaller aphelia distances (see Fig. ^>). These changes de- 
crease the probability of a close encounter between Venus 
and Earth, as can be inferred form Fig.|l|3. Beyond this ke 
value the behaviour of Earth's and Venus' eccentricities do 
not change and in general their orbits turns into a more and 
more regular, quasiperiodic ones as ke increases. Moreover 
their dynamical coupling strengthens. These features can be 
traced on the panels of Fig. 0] 

In the case of Mars the character of the variation of the 
eccentricity is more complex: 

• For the actual masses Mars' eccentricity fluctuates with 
an amplitude of 0.1 and with a very long period of approxi- 
mately 2.5 Myrs. To this fluctuation smaller amplitudes with 
shorter periods are added. It is clear from Fig.QJi, that Mars' 
motion is not coupled to any of the other two terrestrial 
planets. 

• At ke — 2 Mars' eccentricity shows a quite different 
behaviour: it oscillates with a very short period around a 
mean eccentricity of 0.1 with an amplitude of about 0.02. 
With increasing ke the center of oscillation shifts to higher 
values, and reaches its maximum at 5, where the system is 
destabilized. 

• For ke > 5 the motion of Mars gradually turns into a 
regular one. Moreover, dynamical coupling develops between 
all the three inner bodies, which is very well visible from Fig. 

St. 

The terrestrial planets seem to gradually form a subsystem 
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Figure 4. The eccentricities of the inner planets are shown as a function of time for the ke = 1, 2, 15 simulations. The time and mass 
factor evolution of the planets' eccentricities can be traced on the panels. On panel 4c. also Jupiter's and Saturn's eccentricities are 
shown: the existence of the two separate subsystems is well visible. 
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Figure 5. Results of the MEM in the 4 ^ ke ^ 6 region. The 
dotted lines show the ME of the 20 Myr, whereas the solid lines 
show the ME of the 100 Myr integrations. The MEs of Jupiter 
and Saturn are not plotted. The instability window occupies the 
(4.2, 4.9) region. We note that the ordinate is logarithmic. 



as the mass factor increases. As the Earth becomes the dom- 
inant planet in this region and the motions of the two smaller 
neighboring planets are influenced primarily by the Earth 
(the major variation in the eccentricities of the Earth, Venus 
and Mars has the same period, see Fig. [I}:). With regard to 
this we may consider the inner planets as a collection of dy- 
namically mutually dependent planets, namely a subsystem. 

In this mass factor region the giant planets move on 
quiet orbits for the duration of all numerical integrations. 
Jupiter and Saturn also constitute a subsystem which is 
practically not affected from the terrestrial planets. An anal- 
ysis of the eccentricities of all bodies for the subsequent mass 
factors shows that the model graduates into a system which 
consists of two separate and loosely dependent subsystems. 

3.2 The instability window at ke ~ 5 

According to Fig. [3| the ME of Mars reached its maximum 
value at ke — 5. To explore the dynamical evolution of the 
inner planets in more detail simulations were performed in 
the ke — [4,6] region, with a stepsize of 0.1 for 100 Myrs. 
The results of the MEM are shown in Fig. |S] where several 
peaks (ejif = 1-0) in the curve of Mars' ME are visible, which 
correspond to escape orbits. We note that there are two 
minima: at 4.4 and at 4.7. This feature is typical of chaotic 
systems, where small differences in initial conditions, round 
off errors or the applied computing architectures may result 
in different outcomes. To verify this, the two systems with 
ke = 4.4 and 4.7 were integrated on different computers 
and both of the systems decayed. After ke — 5.2 the ME of 
Mars gradually drops down and at ke = 6 it is already less 
than 0.15. 

Fig.|H|shows the evolution of the perihelia (q) and aphe- 
lia (Q) distances of the inner planets for ke ~ 4.6, for which 
case we measured the shortest dynamical lifetime. A strik- 
ing feature is that Mars' perihelion begins to decrease right 
at the start of the integration, as a consequence of the steep 
increase of its eccentricity. At 250 000 yr Mars' e reaches 
0.16 then it decreases to 0.13, around this value it oscillates 
for 1 Myrs. Conversely, q and Q of the Earth and Venus stay 
in a well defined zone. At ~ 1.5 Myr, the perihelion of Mars 
begins to oscillate with quite a large amplitude. The center 
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Figure 6. a.) the inner planets' perihelia (q) and aphelia (Q) 
distances are shown as a function of time for the ke = 4.6 simu- 
lation, b.) The eccentricity of Mars is depicted. Note that Mars' 
orbit is excited to high eccentricity and it becomes Earth-crossing 
at t « 10.5 Myr. c.) The semimajor axes of the Earth and Venus 
are shown together. 



of this oscillation grows secularly and after 10.5 Myr Mars 
becomes an Earth-crosser, moreover it crosses the orbit of 
Venus too. As a consequence of the several close encounters 
with Earth; Mars escapes from the system. This cascade 
mechanism is shown in Fig. Et, where the semimajor axes 
of the Earth and Venus are plotted. 

Because of the chaotic nature of the system, different 
dynamical lifetimes for the different ke values were found. 
Out of the 21 integrations for 100 Myr, Mars escaped four 
times and the system decayed 2 . It must be stressed that the 
chosen length of the integration time is a problem, because 

2 By the decay of the system, we understand that the system has 
evolved into a different one, with a completely distinct configura- 
tion. 
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it may be relatively short for such investigations. It is well- 
known that there exist orbits which are stable for a very 
long time interval and then, all of a sudden they evolve into 
chaotic orbits. These 'sticky orbits' are common in nonlinear 
dynamic al systems and therefore they exist also i n planetary 
systems jDvorak et alll99Stljones fc Sleerl2002>) . This class 
of orbits are embedded in those region of the phase space, 
where regular and irregular orbits are close to each other. 
We suspect that those systems which did not decay, may be 
'sticky systems', and will be unstable at a much longer time. 

According to the chaos theory the extent of this insta- 
bility window is a function of the numerical resolution (in 
our case the Ake). Therefore it is not possible to find the 
exact size of this window, however with our resolution this 
instability window is in the interval 4.2 < ke < 4.9 

The surprising increase in the eccentricity of Mars may 
be a result of some secular resonances between the secular 
frequencies gj and fj. Such phenomena were already found 
by Laskar in the Solar System, who reported that large and 
irregular variations can appear in the eccentricities and in- 
clinations of the terrestrial planets, es pecially of M ercury 
and Mars on time scales of several Gyr jLaskarlll99i) . 

The secular variations of the orbital elements of the 
planets are calculated by means of the Laplace-Lagrange 
theory. Using our MAPLE application the and ifi ' 

amplitudes, the gj and fj secular frequencies and j3j and jj 
were calculated in the ke £ [4,6] interval on a very fine grid 
with a stepsize Akb = 0.001. In this region two frequencies, 
fi and fz have almost the same value (see Table while 
the other ones are well separated. In Table |3] the /s = 0.0 
frequency is not included, and also the corresponding am- 
plitude (2.8402 -10~ 2 ) and the angular phase (106°. 17) were 
left out. 

To visualize the dependence of fi and fz on the mass 
factor, they are shown in Fig. [7Ji, while in Fig. |7J) their 
difference around he = 5 is depicted. From Fig. [7J) it is 
obvious, that the two curves do not intersect each other: the 
difference between the two frequencies are several orders of 
magnitudes greater than the accuracy of the computation, 
which was set to 10~ 10 . The A/ = / 2 — / 3 difference reaches 
its smallest value for ke — 5.03, A/ = 0.1439 " /yr, which 
is at the border of the aforementioned instability window. 

A study of Table [3] shows that the largest amplitude 
are, in the solution for Mars, L 3 2 ^ and for Jupiter 
and and for Saturn , . Accordingly the fi and 
fz frequencies, as was described in section 2.1, can be as- 
sociated with Mars, Jupiter and Saturn. The orbital plane 
of Mars therefore rotates together with those of Jupiter and 
Saturn, giving rise to chaotic behaviour. The equality of two 
apsidal or nodal rates is referred to in Solar System as a sec- 
ular resonance. In this case we have three secular resonances: 
Qm f» fij, fizif ~ f^s and f2j ~ Qs- We suspect that these 
secular resonances are the possible source of the observed 
chaos, and produce the instability window. 



3.3 The 30 < k e ^ 600 region 

Fig. [HJsummarizes the results of the MEM in the 30 ^ ke ^ 
600 region. As one can clearly see from Fig. |S] the MEs of 
the planets gradually increase. None of the curves show any 
peaks, consequently in this region no instability window was 
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Figure 7. On the a.) panel are the fi and / 3 secular frequencies, 
on the b.) panel their difference around 5.0 are plotted versus the 
ke mass factor. 
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Figure 8. Results of the MEM in the 30 < n E < 600 region. 
Note, that the y axis is logarithmic. 

found. These numerical results are in line with the results 
of the linear theory. In this region each of the secular fre- 
quencies have quite different values. Beyond ke ^ 540 Mars 
escaped in all simulations. The MEs of the Earth, Jupiter 
and Saturn behave very similarly and we may say that from 
ke = 100 on the system is dominated by these planets and 
Venus and Mars may be considered as quasiasteroids. 

Generally the semimajor axes of the inner planets are 
confined in a narrow stripe in the order of 10~ 4 -10~ 3 AU. 
These stripes are defined by a very short quasiperiodic os- 
cillation around the mean values of the planets' semimajor 
axes. As the mass factor increases the amplitudes of the os- 
cillations grow, the width of these stripes slightly broadens. 
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Table 3. The fj secular frequencies, the 7^ phases and the amplitudes of the model for ke = 5.03 

calculated using the Lagrange-Laplace theory. The fj are given in arcsec/yr and the ■yjs in degree. 
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Venus 


5.2399 10~ 2 


7.5778 lO" 3 


1.0696 -10~ 2 


2.5795 -10-2 


Earth 


-7.7184 -10- 3 


3.9381 10- 3 


5.9598 -10- 3 


2.3939 -10- 2 


Mars 


2.4032 10~ 3 


-4.2737 lO" 1 


-4.3543 -10- 1 


1.3348 -10- 2 


Jupiter 


4.5852 10~ 6 


3.1885 -10~ 3 


-3.1088 -lO" 3 


-1.2056 -lO" 4 


Saturn 


-2.2595 KT 6 


-7.7642 KT 3 


7.7424 -10~ 3 


-1.8405 -10~ 4 



The typical evolution of the eccentricity of the inner 
planets are shown in Fig. [5] The character of the curves are 
similar to Fig. U, only the periods are shorter. 

Since Venus, Earth and Mars have the same period in 
their eccentricities - just like the Jupiter and Saturn pair 
does - the system is separated into two subsystems: one 
consisting of Venus, Earth and Mars, and the other Jupiter 
and Saturn. This behaviour was already observed for smaller 
mass factor values. The strong coupling of the eccentricities 
of the two giant planets remains unbroken for all Kg. The 
eccentricities of Venus and Mars are mainly determined by 
that of the Earth, which is the superposition of a long period 
(T; ~ 2.5 TO yr) variation with amplitude 0.035 and several 
shorter period (~ 10 4 yr) and smaller amplitude variations 
(see Fig. |5|d) . In the case of Venus, upon this long period 
variation several very short periods are superimposed (see 
Fig. |§Ji) . On the contrary the eccentricity of Mars flickers 
with a period Ti around 0.1 (see Fig.^). 

As the mass factor increases the long period variations 
in the eccentricities of Jupiter and Saturn are more and more 
influenced by the Earth. In order to conserve the total an- 
gular momentum of the system the centers of oscillations of 
Jupiter's and Saturn's eccentricities are anticorrelated with 
the eccentricity of the Earth (see Fig. \10b ) . In spite of the 
large perturbations from the Earth, the coupling between 
the motion of Jupiter and Saturn are still very well deter- 
mined (see Fi e. HOfa). This re s ult is also supported by our 
previous studv lDvorak fc Siilil i2002f) , and therefore we may 
say that the dynamical coupling of Jupiter and Saturn essen- 
tially determines the dynamics of the Solar System's bodies. 

The time developments of the inclinations are very sim- 
ilar to those of the eccentricities: 

• the variations of Jupiter's and Saturn's inclinations are 
coupled, 

• the inclinations of Venus and Mars are primarily deter- 
mined by that of the Earth, 

• with growing ke Jupiter's and Saturn's inclinations are 
increasingly influenced by the Earth. 

When he = 540 Mars' eccentricity had grown secularly 
and at T e ~ 3.4 Myr, Mars passed near by Earth, which 
ejected the planet from the system (see Fig. lllH . In this case, 
we did not observe orbit-crossing. However the minimum 
distance between Earth and Mars was 0.164 AU, which is 
about twice the Hill radius of the massive "Earth" Rh(ke = 
540) = 0. 081 AU. A catacly smic outcome is therefore to be 
expected jjones et al 1 120051) . In the systems with ke = 570 



and 600 Mars was ejected after a sequence of orbit-crossing 
with Earth at T e « 4.35 x 10 6 and at T e « 5.51 x 10 6 years, 
respectively. 

The stability of Venus is an intriguing property of the 
systems. In each of our simulations the orbital elements of 
Venus did not show any sign of chaos. The Earth and Venus 
are prevented from close encounters by a coupling of Venus' 
eccentricity and Azu; Azu is the difference between the lon- 
gitude of perihelion of the Earth and Venus. Whenever the 
perihelion of Earth conjunctions with the aphelion of Venus 
(Aw — 180°), the eccentricity of Venus is around its min- 
imum (see Fig. 1121 . maximizing the distance between the 
two orbits. A similar protection mechanism also occurs with 
real Mars and th e asteroid Pallas. Thi s kind of coupling was 
also observed bv I Jones fc Sleerjl i2002Tl for an 'Earth' in the 
47 UMa system. 



4 DISCUSSION 

We studied the dynamics of a simplified dynamical model 
of the Solar System where we included the three terres- 
trial planets Venus, Earth and Mars and the gas giants 
Jupiter and Saturn. T his model was already studied in detail 
JPvorak fc Sulil2002ft when the masses of all terrestrial plan- 
ets were uniformly increased. It turned out that the different 
systems remained stable for 10 Myrs up to large mass factor 
values. As a continuation of this work the same model with 
a more massive Earth was studied and signs of chaotic be- 
haviour were reported in lDvorak et al.l (I2005T) . In this article 
a detailed exploration of the different S 3 setups is presented, 
and the existence of an instability window is revealed. As an 
important byproduct the plethora of the integrations can 
serve as a general model of exoplanetary systems with two 
massive planets close to the 5:2 mean motion resonance on 
low eccentric orbits for comparable mass ratios of the giant 
planets. 

We investigated over 100 S 3 models where we increased 
the mass of the Earth by a mass factor ke between 1 ^ 
Kg ^ 600. These new systems have been integrated using 
extensive numerical integrations for 20 million years (for se- 
lected systems up to 100 million years) to find out the effect 
of a very massive Earth on the stability of the whole system 
on one hand; on the other hand the model now can serve 
as example of exoplanetary systems of two or three massive 
planets. 

It turned out, that even when the Earth had Jupiter's 
mass (ke ~ 300) and beyond, the system was still stable, but 
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Figure 9. The eccentricities of the Venus, Earth and Mars (from 
top to bottom) for ke = 30 over a 1 million year time-span. 



when ke £ (4.2,4.9), the motion of Mars become chaotic. 
In this instability window Mars' orbital eccentricity finally 
reached values which led to close encounters of Mars with 
the Earth, and even with Venus. After a sequence of close en- 
counters Mars escaped within some millions of years. Using 
the results of the Laplace-Lagrange secular theory we found 
secular resonances acting between the motions of the nodes 
of Mars, Jupiter and Saturn. These secular resonances give 
rise to strong chaos, which is the cause of the appearance of 
the instability window, and eventually the escape of Mars. 

We also found an interesting coupling of Venus' e and 
Azu, which protect Venus form a close encounter with Earth. 
This mechanism was observed in the Solar System and 
in 47 UMa with a hypothetical Earth djones et alJ l2005i 
I Jones fc Sleenll2002l) . 

According to these results, the stability of the Solar Sys- 
tem depends on the masses of the planets, and small changes 
in these parameters may result in a different dynamical evo- 
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Figure 10. a.) The time development of the eccentricities of 
Jupiter and Saturn, b.) The eccentricity of the Earth and the lo- 
cal maximums of Jupiter's and Saturn's eccentricities are shown 
for ke = 250. The coupling between the eccentricities can be well 
observed. 
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Figure 11. This plot shows the escape of Mars for Ke = 540. 



lution of the planetary system. No other instability window 
in ke were found; first results of additional computations 
where we increased the masses of Venus and separately also 
of Mars, showed signs of chaotic motions for some windows 
in k too, but a detailed study of these two other cases of a 
modified S 3 is in preparation. There we also intend to give 
a detailed comparison of the three systems, namely with a 
massive Venus, a massive Earth and a massive Mars. 

Finally we note that all these models may be used as 
dynamical reference models for a better understanding of 
the stability of orbits in extrasolar planetary systems. 



The stability of the terrestrial planets with a more massive "Earth" 11 



0.12 7 
0.10 g 




; - 



90 180 270 360 

Aw [deg] 

Figure 12. This plot shows the relationship between Venus' ec- 
centricity and Avu during 20 Myrs for ke = 540. 
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